Trabajo Final

Trabajo final de la asignatura de Tecnologías para el Análisis de Datos Masivos

true
2021-02-16

1. Contexto

Por cambio climático entendemos tanto el calentamiento global resultado de las emisiones de gases de efecto invernadero, como el incremento en la frecuencia de las catástrofes naturales causadas por las transformaciones en el clima. Este, representa una amenaza existencial y aunque la comunidad científica lo lleva advirtiendo hace años, los gobiernos anteponen políticas con resultados económicamente favorables a corte plazo a determinadas medidas de reducción del cambio climático.

Existen dos problemas destacados por los que no se están aplicando más medidas para combatir al cambio climático:

Con este proyecto se pretende desarrollar una visión intuitiva de la diferencia entre paises en cuanto al uso de energías renovables y convencionales. También como afectan al cambio climático y los efectos negativos este medidos con la diferencia de temperatura.

Los datos tomados se dividen en tres tablas:

La primera tabla se ha obtenido de CDP, Kaggle y se ha contrastado con las fuentes mencionadas para determinar a que hacen referencia exactamente los datos con el fin de no llegar a conclusiones erróneas. La segunda tabla ha sido obtenida directamente de OWID y la tercera tabla de Temperature Change, Kaggle. La tabla de códigos se encuetra en Country Code, Kaggle.

2. Modelo de datos

2.1 Selección de variables

# Carga de datos
co2 = read_csv("Datos/EmisionesC02capita.csv")
energia = read_csv("Datos/ConsumoEnergetico.csv")
temp = read_csv("Datos/CambiosTemperatura.csv")
cod = read_csv("Datos/Codigo_Pais.csv")

Tras la carga de las tablas se han traducido los nombres de las variables.

La tabla de las temperaturas ha requerido de varios cambios para extraer las anomalías de temperatura. En primer lugar se han seleccionado solo las variables útiles para este informe. En segundo lugar, de la variable Month que contiene valores por mes, estación y año, solo se han tomado los valores anuales. La variable Element contiene la desviación estándar y el valor puntual de temperatura por lo que se ha filtrado solo el segundo valor. A continuación, se han pivotado las columnas en las que aparecen los años como variables. Los años se han usado para crear una nueva variable (Año) y los valores de temperatura se han añadido a una columna adicional (Dif_Temp). Esta variable representa la diferencia de temperatura para un año y país concretos con respecto al período de referencia mencionado en le contexto de los datos. Para concluir, se han seleccionadas las variable Pais y Codigo de la tabla de códigos.

Show code
# Selección de variables de cada tabla
co2 %<>% select('Pais' = 'country', 'Codigo' = 'iso_code', 'Año' = 'year',
                'CO2' = 'co2_per_capita')

energia %<>% select('Pais' = 'Entity', 'Codigo' = 'Code', 'Año' = 'Year', 
                    'Carbon' = 'Coal Consumption - EJ', 'Gas' = 'Gas Consumption - EJ', 
                    'Petroleo' = 'Oil Consumption - EJ', 'Nuclear' = 'Nuclear Generation - TWh',
                    'Biomasa' = 'Geo Biomass Other - TWh', 'Hidraulica' = 'Hydro Generation - TWh',
                    'Solar' = 'Solar Generation - TWh', 'Eolica' = 'Wind Generation -TWh')

temp %<>% select('Pais' = 'Area', Months, Element, Y1961:Y2019) %>% 
  filter(Months == 'Meteorological year' & Element == 'Temperature change') %>%
  pivot_longer(c('Y1961':'Y2019'), names_to = 'Año', values_to = 'Dif_Temp') %>% 
  select(Pais, Año, Dif_Temp)

temp$Año = as.numeric(gsub('Y', '', temp$Año))
cod %<>% select('Pais' = 'Country_name', 'Codigo' = 'code_3digit')

Las tres tablas tienen distintos rangos de años, siendo el registro relacionado con la energía el más nuevo por lo que ha marcado el rango de años de las otras dos tablas.

También es diferente el número de países. Con un rápido análisis se ha observado como la tabla de energías solo contiene información de países con una población relevante, mientras que las otras dos monitorizan valores de todo el mundo. Por ejemplo, no se tiene el consumo de energía en la Antártida pero si su temperatura y \(CO_2\). El elevado número de valores de está variable también se debe a que existen valores que agrupan países y a su vez están listados los países que se encuentran en dichas zonas. Es para unificar el nombre de los países, o áreas territoriales, y evitar agrupaciones que se ha añadido la tabla de códigos.

Show code
# Rango de Años y Países
resumen = data.frame(Tabla = c("Energía", "CO2", "Temperatura"),
                     "Inicio_Registro" = c(min(energia$Año), min(co2$Año), min(temp$Año)),
                     "Fin_Registro" = c(max(energia$Año), max(co2$Año), max(temp$Año)),
                     "Países" = c(length(unique(energia$Pais)), length(unique(co2$Pais)), length(unique(temp$Pais)))
                     )
kable(resumen, caption = 'Rango de Años y Total de paises', booktabs = T) %>%
  kable_styling(latex_options = 'hold_position', font_size = 14)
Table 1: Rango de Años y Total de paises
Tabla Inicio_Registro Fin_Registro Países
Energía 1965 2019 100
CO2 1750 2019 236
Temperatura 1961 2019 284

Para pode analizar la relación entre los tres factores (energía, temperatura y \(CO_2\)) en los distintos países, se han unido las tablas mediante un inner join. De esta manera los datos de los países y años resulten estarán completos. En primer lugar, se ha añadido la variable de códigos a las temperaturas y, a continuación, se han añadido las variables de energía y \(CO_2\) mediante la Key: (Código, Año). En segundo lugar, se han seleccionado las variables de interés y, finalmente, se han ordenado las observación manteniendo primero el orden alfabético del nombre del país y luego por orden cronológico.

# Inner Join
datos <- temp %>% inner_join(cod, by='Pais') %>%
  inner_join(energia, by=c('Codigo','Año')) %>% 
  inner_join(co2, by=c('Codigo','Año')) %>%
  select('Pais', 'Codigo', 'Año', 'Dif_Temp','CO2','Carbon':'Eolica')
# Orden cronológico
datos <- datos[order(datos$Año, decreasing=F),] %>% as.data.frame()
Show code
# Tabla Interactiva 1
datos %>% mutate(across(where(is.numeric), ~round(.,3))) %>%
  mutate(across(where(is.character), as.factor)) %>% 
  datatable(rownames = F,
            options = list(pageLength = 5),
            caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;',
                                              'Tabla Interactiva 1: ', 
                                              htmltools::em('Resultado de la unión de las tablas')
   ))

Para concluir la selección de variables se ha cambiado el formato de la información de las energías renovables para contar con una variable categórica. Teniendo en cuenta que cada observación es una combinación de país y año, se han creado las siguientes variables:

A continuación, se han eliminado las variables de generación energética por fuentes renovables.

Show code
# Nuevas variables
renovables = colnames(datos)[10:13]
ren_fun = function(row){
  if (sum(row)!=0) {
    tipo = which.max(row)
    return(renovables[tipo])
  }
  else {
    return("NO")
  }
}

# Renovable
datos$Renovable = apply(datos[,10:13], 1, FUN = ren_fun)

# Gen_Ren
datos %<>% mutate(Gen_Ren = Biomasa + Hidraulica + Solar + Eolica,) %>%
  select(-c("Biomasa", "Hidraulica", "Solar", "Eolica"))

2.2 Limpieza de datos

Valores Nulos

Observando las 5 primeras filas se detecta que hay dos tipos de valores a tener en cuenta, los que son 0 resultantes del registro de información con un valor de nulo, y los NA, resultantes de no registrar esa información específica.

Show code
# Total de NAs
df = data.frame() %>% rbind(rep(0,ncol(datos)))
names(df) = names(datos) 
for (i in 1:ncol(datos)) {
  df[1,i] = sum(is.na(datos[,i]))
}
kable(df, caption = 'Valores no registrados por variable', booktabs = T) %>%
  kable_styling(latex_options = 'hold_position', font_size = 14)
Table 2: Valores no registrados por variable
Pais Codigo Año Dif_Temp CO2 Carbon Gas Petroleo Nuclear Renovable Gen_Ren
0 0 0 207 0 0 0 6 0 0 0

Analizando los países y períodos con NA’s, se ha encontrado una relación entre los períodos de valores faltantes con el resultado de situaciones históricas. Con la disolución en el 1991 de la Unión Soviética y en el 1992 de la República Federal Socialista de Yugoslavia, surgieron varios países como Azerbayan, Estonia o Rusia por un lado y Croacia o Eslovenia de la otra división. Esto explica que haya 12 países con NA’s en la variable de temperatura hasta el 1992. De forma similar Eslovaquia presenta NA’s hasta el 1993 cuando se disolvió Checoslovaquia. Entonces estos hechos históricos tienen un efecto sobre una de las variables, Dif_Temp. Eliminando estas observaciones se perdería la información de las otras variables que si está registrada entonces se han sustituido los NA’s por 0’s, para conservar las otras variables. De ser necesario esto se tendrá en cuenta en futuros análisis.

Otra opción sería analizar los países que hacen frontera con estos problemáticos y hacer una media ponderada de sus valores de temperatura, pero eso requeriría un análisis más exhaustivo teniendo en cuenta la dimensión de cada país y sus respectivas fronteras.

Show code
# Paises con Temperatura no registrada
datos[which(is.na(datos$Dif_Temp)),c('Pais','Año')] -> nans
unique(nans$Pais)
 [1] "Belgium"             "Luxembourg"          "Slovakia"           
 [4] "Singapore"           "Trinidad and Tobago" "Azerbaijan"         
 [7] "Belarus"             "Estonia"             "Kazakhstan"         
[10] "Latvia"              "Lithuania"           "Russia"             
[13] "Turkmenistan"        "Ukraine"             "Uzbekistan"         
[16] "Croatia"             "Slovenia"           

También hay NA’s en el registro de petróleo de Bangladesh hasta el 1971, cuando se independizo de Pakistán. Las variables de energía no representan valores percápita, si no el total del país por lo que no sería correcto asignar los valores de Pakistán a estos fuera de registro. Del mismo modo que han procedido al registrar los datos de las variables de energía de Bangladesh para ese mismo período de tiempo, se han asignado a un valor de 0.

Show code
# Reemplazo de los NA's
datos$Dif_Temp %<>% replace(is.na(datos$Dif_Temp), 0)
datos$Petroleo %<>% replace(is.na(datos$Petroleo), 0)

A raíz del descubrimiento de valores no registrados se ha procedido a comprobar que todos los países tengan el mismo rango de años.

P65 = datos %>% filter(Año == 1965) %>% select(Pais)
P19 = datos %>% filter(Año == 2019) %>% select(Pais)
setdiff(P19$Pais,P65$Pais) 
 [1] "Azerbaijan"   "Belarus"      "Croatia"      "Estonia"     
 [5] "Kazakhstan"   "Latvia"       "Lithuania"    "Russia"      
 [9] "Slovenia"     "Turkmenistan" "Ukraine"      "Uzbekistan"  

Las observaciones de estos países no empiezan en el 1965 como el resto.

Show code
P85 = datos %>% filter(Año == 1985) %>% select(Pais)
URSS <- setdiff(P85$Pais, P65$Pais) #1985 aparecen 10 países de la URSS
P90 = datos %>% filter(Año == 1990) %>% select(Pais)
RFSI <- setdiff(P90$Pais,  P85$Pais) #1990 aparecen dos países de la Rep.Fed.Soc. de Yugoslavia

Se ha confirmado que se empezaron a tener datos de Croacia y Eslovenia en el 1990 y en el 1985 el resto de países.

Valores Atípicos

Se ha considerado que debido a las diferencias de población y territorio entre países, los outliers deben ser detectados por países.

# Detección de Outliers
out_det <- function(data, paises, años, n){ 
  outliers = list()
  for ( pais in paises ) {
    p = data %>% filter(Pais == pais)
    outliers[[pais]] = hampel(p[,n], k = 10, t0 = 3)$ind #Longitud ventana = 2*k+1
  }
  if (isempty(outliers)) {
    return("0 Outliers")
  }
  else{
    df = as.data.frame(plyr::ldply(outliers, rbind))
    colnames(df)[1] = "Pais"
    df %<>%
      pivot_longer(!Pais ,values_to = "indx") %>%
      select(Pais,indx)
  
    for (i in 1:nrow(df)){
      if (df$Pais[i] %in% URSS){
        df$indx[i] = años[df$indx[i + 20]]
      }
      if (df$Pais[i] %in% RFSI){
        df$indx[i] = años[df$indx[i + 25]]
      }
      else {
        df$indx[i] = años[df$indx[i]]
      }
    }
    return(na.omit(df))
  }
}
paises = unique(datos$Pais)
años = unique(datos$Año)

La función out_det() implementa un hample filter para la detección de valores atípicos en series temporales y hace un loop para analizar cada país por separado. Esta se ha aplicado a cada una de la variables.

Show code
# Se ha procedido igual con cada variable
# cada una con su respectiva n, i.e columna
outliers_temp = out_det(datos,paises, años, n = 4)

Todas las variables excepto Gas tienen valores atípicos. A modo de ejemplo se muestran los outliers del consumo de petróleo:

Show code
# Outliers Consumo Petróleo
colnames(outliers_petr) = c("País", "Año")
kable(outliers_petr, caption = 'Outliers Petróleo', booktabs = T) %>%
  kable_styling(full_width = F,latex_options = 'hold_position', font_size = 14) %>%
  column_spec(1,width = "3cm") %>%
  column_spec(2,width = "3cm")
Table 3: Outliers Petróleo
País Año
Finland 1995
Kuwait 1990
Kuwait 1991
Sweden 1979

Todos los valores atípicos han sido analizados y ninguno de ellos es resultado de una medición errónea de las variables. Entonces, se ha decidido no reemplazarlos ya que representan cambios significativos reales que afectan a los países.

La siguiente gráfica muestra la frecuencia de valores atípicos por año.

Show code
# Frecuencia Outliers
outliers <- rbind(outliers_temp, outliers_co2, outliers_carb,
                  outliers_petr, outliers_nuc, outliers_renov)
outliers %>% mutate(Año = as.factor(indx)) %>%
  ggplot(aes(x = Año, fill = Año, color = "tomato1")) +
    geom_bar(stat = 'count', alpha = 0.8) +
    My_Theme + 
    theme(axis.text.x = element_text(angle = 90), legend.position = "none") +
    ggtitle("Frecuencia de Outliers", subtitle = "Por Año")

Como futuros proyectos se podría realizar un análisis con más profundidad de los valores atípicos teniendo en cuenta si son anómalos por encima o por debajo del resto de valores y ver si están relacionados los valores atípicos entre variables.

2.3 Tipología de datos

Show code
# Tipología de datos
datos %<>% mutate(across(Año,as.integer), across(c(Pais,Renovable), as.factor))
str(datos)
'data.frame':   3710 obs. of  11 variables:
 $ Pais     : Factor w/ 72 levels "Algeria","Argentina",..: 1 2 3 4 6 8 9 10 11 12 ...
 $ Codigo   : chr  "DZA" "ARG" "AUS" "AUT" ...
 $ Año      : int  1965 1965 1965 1965 1965 1965 1965 1965 1965 1965 ...
 $ Dif_Temp : num  -0.135 0.063 0.127 -0.882 -0.196 0 -0.079 -0.461 -0.872 0.119 ...
 $ CO2      : num  0.525 2.654 10.683 5.221 0.056 ...
 $ Carbon   : num  0.00293 0.03044 0.7273 0.21281 0 ...
 $ Gas      : num  0.02675 0.148442 0.000106 0.063103 0 ...
 $ Petroleo : num  0.0555 0.9365 0.652 0.2351 0 ...
 $ Nuclear  : num  0 0 0 0 0 0 0 0 0.128 0 ...
 $ Renovable: Factor w/ 5 levels "Biomasa","Eolica",..: 3 3 3 3 4 3 3 3 3 3 ...
 $ Gen_Ren  : num  0.4 1.23 7.99 16.08 0 ...

3. Análisis Exploratorio

En este apartado se ha analizado la distribución de las variables agrupadas por países. A continuación, se ha entrado con más detalle a analizar los valores relacionados con el consumo energético y finalmente se ha analizado la relación de las diferéncias de temperatura con las emisiones de \(CO_2\) percápita.

# Estadísticos
stats <- datos %>% 
  pivot_longer(cols = c(Dif_Temp, CO2, Carbon, Gas, Petroleo, Nuclear, Gen_Ren),
                       names_to = "Variable", values_to = "Valores") %>%
  group_by(Pais, Variable) %>%
  summarise(Media = mean(Valores), Desv.Est. = std(Valores), 
            Min = min(Valores), Q1 = quantile(Valores, 0.25),
            Mediana = median(Valores), Q3 = quantile(Valores, 0.75),
            Max = max(Valores), Sumatorio = sum(Valores))
Show code
# Tabla interactiva 2
stats %>% mutate(across(where(is.numeric), ~round(.,3))) %>%
  mutate(across(where(is.character), as.factor)) %>% 
  datatable(rownames = F, filter = 'bottom',
            options = list(pageLength = 5),
            caption = htmltools::tags$caption(
              style = 'caption-side: bottom; text-align: center;',
              'Tabla Interactiva 2: ',
              htmltools::em('Distribución de las variables contínuas por país, datos del 1965 hasta 2019')
   ))

Por una parte, la media y la desviación típica dan una intuición de como se comporta la función de distribución de las variables. Los valores mínimo, máximo y los cuartiles marcan la distribución del rango de valores de las variables. Y para acabar el sumatorio aporta una medida de la dimensión de la variable durante todo el período del 1965 al 2019. El interés de estos valores en concreto radica en poder comparar los países con la tabla interactiva.

3.1 Consumo energético

Para poder visualizar conjuntamente la energía que proviene de distintas fuentes, se han transforman las energías nuclear y renovables a \(EJ\), donde \(1000TWh = 3.6EJ\).

Show code
# Consumo por país
cons_por_pais19 = datos %>%
  filter(Año == 2019) %>%
  mutate(Nuclear = Nuclear*3.6*10^-3, 
         Renovables =  Gen_Ren*3.6*10^-3,
         `Consumo Total` = Carbon + Gas + Petroleo + Nuclear + Renovables)
# Top 10
top_10 = cons_por_pais19 %>% 
  arrange(-`Consumo Total`) %>%
  head(10) %>% 
  pivot_longer(c( "Consumo Total", 'Carbon':'Nuclear', "Renovables"), 
               names_to = 'Fuente',
               values_to = 'Consumo') %>%
  mutate(Pais = factor(Pais, levels=c(c("China", "United States", "India", "Russia", "Japan", 
                                        "Germany", "Canada", "Saudi Arabia", "Brazil", "Indonesia"))))
top_10 %>%
  ggplot(aes(x = Pais, y = Consumo, fill = Pais)) +
    geom_bar(stat = "identity", alpha = 0.8) +
    facet_wrap(~ Fuente, scales = "free") +
    ylab("Consumo (EJ)") + xlab("País") +
    My_Theme2 + 
   theme(axis.title.x = element_blank(), axis.text.x = element_blank(),axis.ticks.x = element_blank()) +
    ggtitle("Países con mayor consumo energético", subtitle = "Top 10, 2019, por fuente energética primaria")

El conjunto de gráficos muestra como China y EE.UU son los mayores consumidores de energía aunque dependiendo de la fuente primaria hay algunas variaciones en el ranking. El dato más destacable es el caso de Rusia que consume más energía producida con gas natural que China y el hecho que Arabia Saudi se encuentre en el top 10 de consumo total basando su consumo tan solo en petróleo y gas natural.

# Evolución Consumo Total
consumo_total = datos %>%
  group_by(Año) %>% 
  summarise(Carbon = sum(Carbon), Gas = sum(Gas), Petroleo = sum(Petroleo),
            Nuclear = sum(Nuclear)*3.6*10^-3, Renovables = sum(Gen_Ren)*3.6*10^-3)
Show code
consumo_total %>%
  pivot_longer(c('Carbon':'Renovables'), 
               names_to = 'Fuente',
               values_to = 'Consumo') %>% 
  mutate(Fuente = factor(Fuente, levels=c("Renovables","Nuclear","Gas","Petroleo","Carbon"))) %>%
    ggplot(aes(x = Año, y = Consumo, fill = Fuente)) +
      geom_area(alpha = 0.8) + 
      scale_x_continuous(breaks = seq(1965,2020,5)) +
      ylab("Consumo(EJ)") +
      ggtitle("Consumo Energético Total(EJ), 1965 a 2019") + 
      My_Theme -> plot3
ggplotly(plot3)

En este gráfico se puede ver como el consumo energético total se ha más que cuadriplicado en las últimas 6 décadas, especificamente para los países contenidos en estos datos, de 118.671 a 477.48 Exajoules. Este incremento coincide con el aumento de la población mundial que durante este período, se ha duplicado, y con que la mayoría de países han aumentado sus niveles de consumo.

# Evolución Consumo proporcional
consumo_prop = consumo_total %>%
  mutate(Año, sum = Carbon + Gas + Petroleo + Nuclear + Renovables,
                        Carbon = 100*Carbon/sum,Gas = 100*Gas/sum,
                        Petroleo = 100*Petroleo/sum,Nuclear = 100*Nuclear/sum, 
                        Renovables = 100*Renovables/sum) %>% select(-c(sum))

El siguiente gráfico representa la evolución de la proporción de energía consumida de todos los países, según la fuente primaria.

Show code
consumo_prop %>%
  pivot_longer(c('Carbon':'Renovables'), names_to = 'Fuente', values_to = 'Consumo') %>%
  mutate(Fuente =factor(Fuente, levels=c("Renovables","Nuclear","Gas","Petroleo","Carbon"))) %>%
  ggplot(aes(x = Año, y = Consumo)) +
    geom_line(aes(color=Fuente),size = 1.5, alpha=0.7) +
    scale_x_continuous(breaks=seq(1965,2020,5)) +
    ylab("Consumo%") +
    ggtitle("Proporción de consumo energético por fuente primaria",
            subtitle = "Entre los años 1965 y 2019")+
    My_Theme -> plot4
ggplotly(plot4)

En general el uso del petróleo ha disminuido en las últimas 3 décadas. El uso de gas natural como fuente de energía ha aumentado y el del carbón ha disminuido con respecto a los años 60. La energía nuclear y la que proveniene de fuentes renovables siguen siendo las más bajas aunque han aumentado un 6.66% y un 153.1% respectivamente en los últimos 20 años.

3.2 Temperatura y CO2

Varios estudios muestran la relación entre las emisiones de \(CO_2\) y el aumento de la temperatura global de la tierra. Pero, ¿hay relación entre las emisiones de \(CO_2\) percápita de un país y los cambios de temperatura que sufre?

Los siguientes valores son el promedio entre todos los países del conjunto de datos.

Show code
datos %>% group_by(Año) %>%
  summarise(meanco2 = round(mean(CO2),2), meantemp = round(mean(Dif_Temp),2)) %>%
  ggplot(aes(y=meantemp, x=Año)) +
    geom_point(aes(color=meanco2), alpha = 0.8, size=3.5) +
    geom_line(color="#000099") + 
    geom_hline(yintercept = 0, linetype = "dashed", color = "#FF0000") +
    annotate(geom = "text", y = -0.2, x = 2010, label = "Respecto el período del 1951 al 1980", color="#FF0000") + 
    scale_colour_gradientn(colours=rainbow(4)) + 
    scale_x_continuous(breaks=seq(1965,2020,5)) +
    labs(title = "Evolución de la Diferencia de Temperatura promedio",
        y = "Temperatura (ºC)",
        color = "Emisiones per capita (kg)") +
    My_Theme -> emisiones
ggplotly(emisiones)  

Se puede observar como hasta mediados de los 80 la temperatura era más o menos estable. A partir de los 90, a excepción de algunos años esta temperatura no ha hecho más que aumentar.

# CO2 Vs Diferencias de temperatura
datos %>% group_by(Pais) %>% 
  summarise(CO2 = round(mean(CO2),3), Dif_Temp = round(mean(Dif_Temp),3)) %>%
  ggplot(aes(y=CO2, x=Dif_Temp, color = Pais)) +
    geom_point(size= 4, alpha = 0.8) +
    theme(legend.position= "none") +
    labs(title="CO2 Vs Diferencia Temperatura") +
    My_Theme -> co2_vs_temp

ggplotly(co2_vs_temp)

No existe ninguna relación directa entre las anomalías de temperatura y las emisiones de \(CO_2\) por persona. Analizando la distribución de los países respecto al eje horizontal parece que la diferencia media de temperaturas está relacionado con la ubicación geográfica de los países. En el siguiente enlace hay un gráfico que muestra la evolución del cambio de de temperaturas en la superficie terrestre sobre un mapa mundial y se observar como se desplaza en forma de varías franjas de temperatura que cubre varios países simultaneamente.

4. Aprendizaje Estadístico

En este apartado se han aplicado métodos de aprendizaje estadístico a los datos analizados previamente. Inicialmente se ha realizado un análisis de componentes principales que junto con el clustering han representado dos formas distintas de agrupar los países. En segundo lugar, se han comparado el algoritmo SVM y el Random Forest para tratar de clasificar los países según la fuente de energía renovable con la que más energía generan. Finalmente, se ha comparado la eficiencia de una RLM implementando un método selección de variables (Forward Stepwise) con la de un Random Forest, esta vez en la predicción de las variaciones de temperatura de la superficie terrestre.

4.1 Analisis de Componentes Principales

Para el siguiente análisis se han escogido los datos del año más reciente, 2019, para ser más próximo a la realidad actual.

# ACP
datos %>% filter(Año == max(Año)) %>%
  select(-c(Pais, Codigo, Año)) -> datos_acp
paises = datos %>% filter(Año == max(Año)) %>% select(Codigo)
rownames(datos_acp) = paises$Codigo
# ONE-HOT Encoding
dummy = dummyVars(~ Renovable, data = datos_acp )
Renovable_dummy = as.data.frame(predict(dummy, newdata = datos_acp))
Renovable_dummy %<>% select(c(Biomasa = Renovable.Biomasa, Eolica = Renovable.Eolica,
                             Hidraulica = Renovable.Hidraulica, Solar = Renovable.Solar))
# Ensamblaje
datos_acp %<>% cbind(Renovable_dummy) %>% select(-c(Renovable))

cp <- prcomp(datos_acp, scale = T)
Show code
#Variación Explicada Acumulada
cp.var = as.data.frame(cp$sdev^2) #Varianzas
cp.expl = cp.var/sum(cp.var) #proporción de variación explicada
colnames(cp.expl) = "var_explicada"
ggplot(data = cumsum(cp.expl), aes(var_explicada, x=1:11))+
  geom_line(color='green2') + geom_point() +
  geom_hline(yintercept=1, linetype='dashed', color= "#FF3333") +
  geom_bar(stat = "identity", alpha = 0.8, fill = "#99FFFF")+
  scale_x_continuous(breaks=seq(1,11,1)) +
  ggtitle('Variación Explicada Acumulada') +
  xlab('Componentes Principales') + ylab('Proporción de Variación') +
  My_Theme2

Se puede afirmar entonces que las 4 primeras componentes principales explican más de un \(75\%\) de la variabilidad de los datos. Respecto estas componentes, en la tabla 4 se pueden ver sus cargas que describen la relación entre las respectivas componentes y las variables originales.

Las otras componentes ya son una mezcla en la se observa parte de la contribución de la temperatura y una alta correlación con tecnologías renovables menos frecuentes que la hidráulica.

Show code
# Cargas
cargas = as.data.frame(cp$rotation[,1:4])

kableExtra::kable(t(cargas), caption = 'Cargas de los 4 primeros componentes principales', booktabs = T) %>%
  kable_styling(latex_options = 'hold_position', font_size = 20)
Table 4: Cargas de los 4 primeros componentes principales
Dif_Temp CO2 Carbon Gas Petroleo Nuclear Gen_Ren Biomasa Eolica Hidraulica Solar
PC1 0.0844983 -0.0522465 -0.3882330 -0.4373038 -0.4948476 -0.4403145 -0.4481719 0.0654685 -0.0580734 -0.0094363 0.0285337
PC2 -0.0534866 -0.5014625 0.1294046 -0.1129644 -0.0487856 -0.0553438 0.1202391 -0.2283824 -0.2877482 0.6381544 -0.3925163
PC3 0.3615064 -0.2743393 -0.1034509 0.0273478 -0.0061710 0.1044472 -0.0618161 -0.0806431 0.6808475 -0.1783050 -0.5157356
PC4 -0.2526690 -0.0988651 -0.1608380 0.0667643 0.0008057 0.0250724 -0.1160086 -0.8199820 0.2513624 0.0244250 0.3818643

En el siguiente gráfico se puede evaluar la estructura de los países. Respecto a la primera componente destacan dos grandes aglomeraciones a la derecha y dos puntos atípicos con valores muy negativos. Tras el análisis de las CP1 y CP2 se entienden los valores de China y EE.UU debido al alto consumo energético de ambas. Luego la separación respecto la segunda componente es menos directa ya que por ejemplo separa Italia de España que producen unas emisiones de \(CO_2/capita\) similares. El hecho por el que se han representado divididas es por la contribución del otro factor de la CP2, es decir Italia emplea la energía hidráulica como su principal fuente de energía renovable, mientras que en España predomina la eólica. Este hecho, explica la mayor parte de la división entre los dos segmentos cerca del eje de ordenadas.

Puede ser interesante realizar consultas con la tabla interactiva 2 para entender mejor esta segmentación de los países.

Show code
# Proyección sobre CP1-CP2
proyec = as.data.frame(cp$x[,1:2])
proyec %<>% setDT(keep.rownames=T)

ggplot(proyec, aes(x = PC1, y = PC2))+
  geom_point(size = 1.8, color = "tomato1", alpha = 0.8) +
  geom_text_repel(aes(label = rn), size = 2, max.overlaps = 30) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") + 
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
  xlim(-11,4) + 
  ggtitle("Proyección de los Paises sobre la CP1 y CP2") +
  My_Theme2

Show code
#Filtro Códigos
datos %>% filter(Año == max(Año)) %>%
  select(Pais, Codigo) %>% 
  mutate(Codigo = as.factor(Codigo)) %>%
  datatable(rownames = F, filter = 'top',
            options = list(pageLength = 1, dom = "t"),
            caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;',
                                              'Tabla Interactiva 3', 
                                              htmltools::em('Filtra el país del cual quieras saber su código o viceversa'))
            )

Sería interesante para proyectos centrados en ACP, comparar la evolución de las componentes principales, i.e. como cambia el foco de variabilidad entre países con los año. También añadir otras variables para obtener indices que describan el desarrollo de los países o el nivel de adversidad frente a los efecto del cambio climático.

4.2 Clustering de países

K-Means

Siguiendo con la tónica del análisis anterior se ha desarrollado una segmentación de los países usando el algoritmo KMeans.

# KMeans
datos_num = r.datos_acp #datos creados para la acp, con los dummys

from sklearn.cluster import KMeans
sdc = [] #suma de distancias cuadradas al centro del cluster respectivo
for i in range(1,12):
  kmeans = KMeans(n_clusters = i, init = "k-means++", max_iter = 300, n_init = 10, random_state = 0)
  kmeans.fit(datos_num)
  sdc.append(kmeans.inertia_)
Show code
# Visualización del error por k
_ = plt.figure() 
_ = plt.plot(range(1,12), sdc, "o--")
_ = plt.title("Suma de distáncias cuadradas")
_ = plt.xlabel("Número de Clusters(k)")
_ = plt.ylabel("Suma")
plt.show()

# Ajustar los datos con KMeans
kmeans = KMeans(n_clusters = 4, init="k-means++", max_iter = 300, n_init = 10, random_state = 0)
y_kmeans = kmeans.fit_predict(datos_num)

Para aproximarlo a la representación del apartado anterior, se han usado como dimensiones la variable Petroleo(CP1) y la \(CO_2\)(CP2).

Show code
# Visualización de la segmentación
_ = plt.figure() 
labels = ["Bajo Consumo de Petróleo", "CHN", "USA", "Consumo Moderado de Petróleo"]
colors = ["red", "blue", "green", "cyan"]
for i in range(0,4):
  x = datos_num.Petroleo[y_kmeans == i]
  y = datos_num.CO2[y_kmeans == i]
  if len(x) > 50:
    _ = plt.scatter(datos_num.Petroleo[y_kmeans == i], datos_num.CO2[y_kmeans == i], s = 100, c = colors[i], label = labels[i], alpha=0.2)
  else: 
    _ = plt.scatter(datos_num.Petroleo[y_kmeans == i], datos_num.CO2[y_kmeans == i], s = 100, c = colors[i], label = labels[i])
_ = plt.scatter(kmeans.cluster_centers_[:,4], kmeans.cluster_centers_[:,1], s = 10, c = "yellow", label = "Baricentros")
_ = plt.title("Cluster de países")
_ = plt.xlabel("Consumo Energía fuente Petroleo(EJ)")
_ = plt.ylabel("Emisiones CO2 (kg/capita)")
_ = plt.legend()
plt.show()

Con la regla del codo se ha seleccionado \(k=4\) como el número óptimo de clusters entre los países de estudio. Este valor coincidiendo con la segmentación del ACP. Como en el apartado anterior, se observan China y EE.UU en sus propias burbujas lejanos de cualquier otro país, hecho que no es de extrañar y probablemente aplicable en muchas otras métricas. Con este modelo no se ve tan clara la división de las dos aglomeraciones. A grandes rasgosc el modelo ha vuelto a dividir por la variable Petroleo, y algún otro factor no representado en el gráfico como por ejemplo la variable Hidraulica que con el ACP se ha visto que es un factor importante en la segmentación de los países.

4.3 Clasificación de la fuente de renovable dominante

Los datos de estudio contienen dos variables categóricas, Pais y Renovable. Debido al elevado número de valores de la primera variable, 72 países distintos, se ha seleccionado la segunda con 5 valores como target de la clasificación.

Para entrenar este modelo no se han usado las variables de generación energética de las fuentes renovables. Esto se debe a que la variable dependiente de esta clasificación es función de estas, concretamente el máximo de estas.

Destacar que para este apartado y el siguiente, al tratarse de algoritmos de aprendizaje supervisado se ha usado todo el conjunto de daots, seperando una parte de este para el entrenamiento y otra para el test de los modelo.

Show code
# Valores de la variable objetivo 65-19
datos_clas <- datos %>% select(c(Año, Dif_Temp, CO2, Carbon, Petroleo, Gas, Nuclear, Gen_Ren, Renovable))
frec_ren <- as.data.frame(prop.table(table(datos_clas$Renovable)))
df = data.frame() %>% rbind(round(frec_ren$Freq*100,2))
colnames(df) = frec_ren$Var1
kable(df, caption = 'Distribución de valores(%), 1965-2019') %>%
  kable_styling(latex_options = 'hold_position', font_size = 14)
Table 5: Distribución de valores(%), 1965-2019
Biomasa Eolica Hidraulica NO Solar
6.74 4.1 78.14 9.11 1.91

Para equilibrar un poco las clases se han cogido datos a partir del 2013 que es cuando ya todos los países usan alguna fuente de energía renovable, i.e. desaparece la categoría NO, y las observaciones se distribuyen un poco más en las distintas categorías, como muestra la Tabla 6.

Show code
# Valores de la variable objetivo 13-19
datos_clas %<>% filter(Año>=2013)
frec_ren <- as.data.frame(prop.table(table(datos_clas$Renovable)))
df = data.frame() %>% rbind(round(frec_ren$Freq*100,2))
colnames(df) = frec_ren$Var1
kable(df, caption = 'Distribución de valores(%), 2013-2019') %>%
  kable_styling(latex_options = 'hold_position', font_size = 14)
Table 6: Distribución de valores(%), 2013-2019
Biomasa Eolica Hidraulica NO Solar
9.92 17.26 62.5 0 10.32

Este es un problema de multiclasificación con varias dimensiones. Se han usado dos algoritmos distinos, el SVM, support vector machine y el Random Forest. El siguiente código muestra el preprocesamiento de los datos para ambos modelos.

# Separación de variables
X = r.datos_clas.iloc[:, 0:-1]
y = r.datos_clas.iloc[:,-1]

#Train/Test split (75/25)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state = 0)

#Estandarización 
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

Support Vector Machine

Se han entrenado dos modelos de SVM uno con kernel rbf y otro polinomial de grado 3. Seguidamente se ha medido su eficiencia clasificando el set de entrenamiento mediante dos métricas, el porcentaje de observaciones correctamente clasificadas(accuracy) y el valor f1.

# Ajustar el SVM en el Conjunto de Entrenamiento
from sklearn.svm import SVC
rbf = SVC(kernel = 'rbf', random_state = 1, gamma = "scale")
rbf_pred = rbf.fit(X_train, y_train).predict(X_test)
poly = SVC(kernel = "poly", random_state = 1, gamma = "scale")
poly_pred = poly.fit(X_train, y_train).predict(X_test)

# Evaluación de los modelos
from sklearn.metrics import accuracy_score
from sklearn.metrics import f1_score
from sklearn.metrics import confusion_matrix
rbf_accuracy = accuracy_score(y_test, rbf_pred)
poly_accuracy = accuracy_score(y_test, poly_pred)
# rbf_f1 = f1_score(y_test, rbf_pred, average='weighted')
# poly_f1 = f1_score(y_test, poly_pred, average='weighted')

Las dos ultimas lineas de código están comentadas, pero se han ejecutado en una chunk escondida debido a un aviso inevitable de Python de la función f1_score.

Show code
# Precisión
print("", 'Accuracy (RBF Kernel): ', "%.2f" % (rbf_accuracy*100), '\n',
  'F1 (RBF Kernel): ', "%.2f" % (rbf_f1*100), '\n',
  'Accuracy (Poly Kernel): ', "%.2f" % (poly_accuracy*100),'\n',
  'F1 (Poly Kernel): ', "%.2f" % (poly_f1*100))
 Accuracy (RBF Kernel):  72.22 
 F1 (RBF Kernel):  63.41 
 Accuracy (Poly Kernel):  72.22 
 F1 (Poly Kernel):  62.90
Show code
# Matriz de Confusión SVM
cm = confusion_matrix(y_test, poly_pred)
df_cm = pd.DataFrame(cm, index = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]],
                  columns = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]])
_ = plt.figure(figsize = (10,7))
_ = sn.heatmap(df_cm, annot=True)
_ = plt.title("Matriz de Confusión")
_ = plt.xlabel("Predicción")
_ = plt.ylabel("Real")
plt.show()

Random Forest

Se ha entrenado un modelo y se ha evaluado de la misma manera que el SVM.

# Ajustar un RF a los datos de train
from sklearn.ensemble import RandomForestClassifier
rf_model = RandomForestClassifier(n_estimators = 500,
                                  criterion = "entropy", 
                                  random_state = 0) 
# Predicción
rf_pred = rf_model.fit(X_train, y_train).predict(X_test)

# Evaluación
rf_accuracy = accuracy_score(y_test,rf_pred)
rf_f1 = f1_score(y_test, rf_pred, average='weighted')
Show code
# Precisión
print("", 'Accuracy (RBF Kernel): ', "%.2f" % (rf_accuracy*100),'\n',
      'F1 (RBF Kernel): ', "%.2f" % (rf_f1*100))
 Accuracy (RBF Kernel):  92.86 
 F1 (RBF Kernel):  92.95
Show code
# Matriz de Confusion RF
cm = confusion_matrix(y_test, rf_pred)
df_cm = pd.DataFrame(cm, index = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]],
                  columns = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]])
_ = plt.figure(figsize = (10,7))
_ = sn.heatmap(df_cm, annot=True)
_ = plt.title("Matriz de Confusión Random Forest")
_ = plt.xlabel("Predicción")
_ = plt.ylabel("Real")
plt.show()

En la precisión de los dos modelos vemos que la del Random forest es significativamente mejor. Esto se ve representado en las matrices de confusión, donde vemos que las dos predicen igual la clase Hidráulica pero el segundo modelo aumenta la precisión en las otras tres clases.

4.4 Predicción de las anomalías de Temperatura

Para poder aplicar algoritmos de predicción a los datos de estudio se ha decidido predecir los valores de Dif_Temp apartir de las otras variables. En la realidad esto no tendría mucho sentido ya que las variables son registradas de forma simultanea por lo que no tiene valor predecir la anomalía de temperatura que ocurre en el mismo momento en que se puede medir.

RLM con selección Forwars Stepwise

# Selección de Variables
datos_reg = datos %>% select(-c(Pais, Codigo))

#ONE-HOT Encoding
dummy_reg = dummyVars(~ Renovable, data = datos_reg)
Renovable_dummy = as.data.frame(predict(dummy_reg, newdata = datos_reg))
Renovable_dummy %<>% select(c(Biomasa = Renovable.Biomasa, Eolica = Renovable.Eolica,
                             Hidraulica = Renovable.Hidraulica, Solar = Renovable.Solar))
#Ensamblaje
datos_reg %<>% cbind(Renovable_dummy) %>% select(-c(Renovable))

set.seed(12)

#Train/Test 75/25%
train_idx = sample(nrow(datos_reg), nrow(datos_reg)*0.75)
datos.train = datos_reg[train_idx, ]
datos.test = datos_reg[-train_idx, ]

num_var <- as.numeric(ncol(datos_acp) - 1) #Vars explicativas

# Función de predicción
predict = function(object, newdata, id, ...){
  form = as.formula(object$call[[2]])
  mat = model.matrix(form, newdata)
  coefi = coef(object,id=id)
  xvars = names(coefi)
  mat[,xvars]%*%coefi
}

# 10folds CV
k = 10
folds = sample(1:k, nrow(datos.train), replace=TRUE)
cv.errors = matrix(NA,k,num_var, dimnames =list(NULL , paste(1:num_var)))
for(j in 1:k){
  best.fit = regsubsets(Dif_Temp~.,data = datos.train[folds!=j,],
                      nvmax = num_var, method='forward') #Forward SW
  for(i in 1:num_var){
    pred = predict(best.fit, datos.train[folds==j,],id = i)
    cv.errors[j,i]=mean( (datos.train$Dif_Temp[folds==j]-pred)^2)
  }
}
rmse.cv = sqrt(apply(cv.errors, 2, mean))
plot(rmse.cv, pch=19, type="b", main="Error", xlab="Número de variables del modelo")
# Conjunto óptimo de variables
reg.sub = regsubsets(Dif_Temp~.,data = datos.train, nvmax = num_var, method='forward')
coef(reg.sub, id = 4)
  (Intercept)           Año      Petroleo       Nuclear        Eolica 
-61.978279733   0.031379016  -0.012744912   0.000667009   0.210211552 

Estas son las 4 variables que forman el modelo óptimo encontrado siguiendo el método de selección forward step-wise. Seguidamente se ha calculado una regresión lineal múltiple con estas variables.

# Modelo Final
reg.f = lm(Dif_Temp ~ Año + Petroleo + Nuclear + Eolica, data = datos.train)
lm.pred = predict(reg.f, newdata = datos.test)
error.rlm = mean((datos.test[, "Dif_Temp"] - lm.pred)^2) 
summary(reg.f)

Call:
lm(formula = Dif_Temp ~ Año + Petroleo + Nuclear + Eolica, data = datos.train)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.02171 -0.34399 -0.01749  0.32684  1.91119 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.198e+01  1.312e+00 -47.230  < 2e-16 ***
Año          3.138e-02  6.588e-04  47.633  < 2e-16 ***
Petroleo    -1.274e-02  3.341e-03  -3.815 0.000139 ***
Nuclear      6.670e-04  1.748e-04   3.815 0.000139 ***
Eolica       2.102e-01  5.354e-02   3.926 8.84e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5285 on 2777 degrees of freedom
Multiple R-squared:  0.4851,    Adjusted R-squared:  0.4843 
F-statistic:   654 on 4 and 2777 DF,  p-value: < 2.2e-16

De la regresión anterior todas las variables son estadisticamente significativas al \(1\%\). Respecto la precisión del modelo para explicar las diferencia de temperaturas, el \(R^2_{ajustado}\) indica que la regresión explica alrededor de la mitad de variable dependiente.

Random Forest Regression

Ya que este algoritmo ha sido tan eficiente en la clasificación del apartado anterior, ha sido seleccionado también en la predicción de las diferencias de temperatura. Para entrenar el Random Forest se van han usado las mismas variables y conjuntos de entrenamiento y test que con la RLM.

Mediante un grid searching se ha calculado el valor optimo para el parámetro m del modelo. Este parámetro representa el número de variables, elegidas aleatoriamente, usadas como candidatos en cada split).

Show code
# Tuning de m
mtry = c(2,4,6,8,10)
errors = matrix(NA , length(mtry) , 1, dimnames =list(paste(mtry), NULL))
set.seed(12)
for(i in 1:length(mtry)){
  rf = randomForest(Dif_Temp~ . , data = datos.train, mtry = mtry[i],
                 metric = 'RMSE')
  errors[i,1] = mean(rf$rsq)
}
plot(x = rownames(errors), errors, pch=19, type = "b", 
     main = "Error", xlab= "m", ylab= "RMSE")

Show code
# Modelo Final
rf.f = randomForest(Dif_Temp ~ . ,
                    data = datos.train,
                    ntree = 500, mtry = 2,
                    metric = 'RMSE', importance = T)
# rf.pred = predict(rf.f, newdata = datos.test)
# error.rf = mean((datos.test[, "Dif_Temp"] - rf.pred)^2)

#Importancias de las variables 
rf.f$importance %>% as.data.frame() %>%
  select(Prop_ECM = "%IncMSE") %>%
  cbind(Var = rownames(rf.f$importance)) %>%
  ggplot(aes(x = Var, y = Prop_ECM)) +
    geom_bar(stat = "identity", aes(fill = Var)) +
    theme_minimal() +
    theme(axis.title.x = element_blank(), axis.text.x = element_blank(),axis.ticks.x = element_blank()) +
    ggtitle('Proporción de reducción Media del ECM')

Un estudio interesante sobre la predicción de anomalías de temperatura que se podría hacer con estos datos es una predicción por serie temporal. En este caso si se estarían prediciendo valores que aún serían desconocidos.

Conclusiones

El objetivo del proyecto era desarrollar una comparativa entre el consumo energético de diferentes países y sus emisiones de \(CO_2\), además de tener una visión de las anomalías de temperatura que estamos sufriendo en los últimos años a nivel global. Creo que con el anterior análisis se ha alcanzado este objetivos. Además, han aparecido nuevas idea más concretas sobre que caminos de investigación pueden ser interesantes en este ámbito. Uno que me interesa en especial es el desarrollo de indicadores para comparar los países. Otro ejemplo sería el análisis de la relación entre las emisiones de \(CO_2\) percápita y la fuente de energía primaria más usada en el país. En conclusión, a mi parecer este ha sido un proyecto muy interesante con el que he aprendido mucho no solo en conceptos, también en intereses profesionales.